This is a model to use AoT network data for calculating the ambient temperatures in deifferent parts of an urban area. It uses data from AoT sensors and from Land Surface Temperature calculated from remote sensing model to train a model for calculatting ambient temperature. Then it uses that temperature with calculated ambient humidity from spatial statistics to calculate heatwave index in each part (probably a 30m x 30m area) of the area and raise a flag for areas experiencing higher than threshold. It can later be used for the following two purposes: * Adaptation Planning: a predictive model that can show potential areas of high heatwave index. * Mitigation Planning: a seperate classification model for suggesting the best heatwave mitigation Urban Planning strategy for each specific place in the area.

Required (available) data sets:

Start by looking into the location of nodes throughout Chicago

library(sp) #spatial data wrangling & analysis

library(rgdal) #spatial data wrangling
library(rgeos) #spatial data wrangling & analytics
library(tidyverse) # data wrangling

#Getting the data of nodes
load("nodes.rda")

Now Getting the total population in each Census Block Group based on ACS-2016 5-year estimate

library(tidycensus)
library(tidyverse)
library(sf)

load("ACS_population.rda")

Mapping Nodes on Block Groups

Here there is a basic map showing the locations of nodes and the Census Block Groups

library(mapview)
library(RColorBrewer)

# Blocks
m1 <- mapview(chicago)

# nodes
m2 <- mapview(nodes.spt, size=5, legend=FALSE, color="red", popup=popupTable(nodes.spt, zcol = c("address", "node_id")))

m1+m2

It’s look like the avialable nodes do not cover all the Census Block. To reduce the area of study, I’ll add Chicago community areas as a layer.

# First to remove the outlier node!
outlier.index <- which.min(nodes.spt$lon)
nodes.spt <- nodes.spt[-outlier.index,]
Stations <- Stations[-outlier.index,]
temp.daily <- temp.daily[,,-outlier.index]


# Now to eliminate BGs outside the containing box
this.area <-st_sfc(st_polygon(list(rbind(c(min(nodes.spt$lon), min(nodes.spt$lat)),
                             c(max(nodes.spt$lon), min(nodes.spt$lat)),
                             c(max(nodes.spt$lon), max(nodes.spt$lat)),
                             c(min(nodes.spt$lon), max(nodes.spt$lat)),
                             c(min(nodes.spt$lon), min(nodes.spt$lat)))))) %>% st_sf(relation_to_geometry = NA_character_)

st_crs(this.area) <- st_crs(chicago)
#st_crs(this.area) <- "+proj=longlat +datum=NAD83 +no_defs"
new_chicago <- st_intersection(this.area, chicago)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
m2 <- mapview(nodes.spt, size=5, legend=FALSE, color="red", popup=popupTable(nodes.spt, zcol = c("address", "node_id")))

m3 <- mapview(new_chicago)
m2 + m3

## Getting LandSat 8 data and calculate Land Surface Temperature (LST) The Landsat is a multispectral satellite currently on the eighth of it’s series. Landsat-8 data is freely available on the USGS’s Earth Explorer website. All we need to do is sign up and find a scene that match our area of study. The notation used to catalog Landsat-8 images is called Worldwide Reference System 2 (WRS-2). The Landsat follows the same paths imaging the earth every 16 days. Each path is split into multiple rows. So, each scene has a path and a row. 16 days later, another scene will have the same path and row than the previous scene. This is the essence of the WRS-2 system.

USGS provides Shape files of these paths and rows that let us quickly visualize, interact and select the important images. After checking the shape file of Landsat 4-8, WRS-2, I found that city of Chicago is within Path 23 and row 31. ( or can use a converter to find Path and Row from latitude and longitude).

at this time I used GloVis tool to download 6 images of Landsat 8 OLI/TIRS C1 Level-1. Then to test the model and for the proof of concept I am just adding one layer (band 6) of downloaded LandSat GeoTiff Image.

source("LandSat_band6.R")
mapviewOptions(mapview.maxpixels = 5e+05)
m4 <- mapview(landsat, legend=FALSE)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 62343791 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  62343791 '
m4+m2+m3

We need to crop the raster part to match the area of study:

# landsat <- crop(landsat, new_chicago)